/*
 *    /||\
 *   / || \
 *  /__||__\
 *   __
 *  \  | UNIVERSIDADE DE PASSO FUNDO
 *   \ | Instituto de Ciências Exatas e Geociências
 *    \| Curso de Ciência da Computação
 * 
 * Trabalho de Conclusão II - 2008/2
 * 
 * Título do trabalho: Paralelização do Solver Intervalar LSS
 * Acadêmico: Alexandre Almeida
 * Orientador: Prof. Dr. Carlos Amaral Hölbig
 * 
 * Módulo: CalculateDR2temp.h / CalculateDR2temp.cpp
 * Descrição: Implementa a versão paralela do cálculo das operações 3 e 4 da Tabela xxx do TCC.
 * 
 */

#include "CalculateDR2temp.h"

/**
 * calculateDR2temp - Function called by the MasterLss class. This function computes (Ai^-1 = I^-1 * A^-1) and (Ea = (I^-1 * A^-1) - Ai^-1),
 *                    where I^-1 = R2 parameter, and A^-1 = R1 parameter. This function distributes the rmatrix "R2" among the processors in a
 *                    stripped partitioning fashion, the same way IminusAB function does, and gathers the R1 matrix in every processor, using an MPI_Allgather operation.
 *                    This Allgather operation makes use of the already distributed R1 matrix, that was split among the processors in the IminusAB
 *                    function.
 *                    Also, each slave, and the master itself, can compute its partial results using threads, depending on the parameter "numThreads".
 *                    If it's greater than 1, then the program will subdivede the domain into "numThreads" threads.
 *            
 * Input parameters:
 *      rmatrix& R2: left operand of A * B;
 *      rmatrix& R1: right operand of A * B.
 * 
 * Output parameters:
 *       rmatrix& D: result of A * B;
 *  rmatrix& R2temp: result of (A * B) - D, where (A * B) is computed using the long accumulator.
 * 
 * Return value: void.
 * 
 */
void calculateDR2temp(rmatrix& R2, rmatrix& R1, rmatrix& R2temp, rmatrix& D, int rank, int numThreads, int worldCommSize){
	// Calculates the total size of R1, so the master can broadcast it to the slaves,
	// each of which will gather a copy of R1.
	int sizeR1 = Ub(R1, ROW) - Lb(R1, ROW) + 1;
	// Determines the part size of R2 that will be sent to each process
	int slice = sizeR1 / worldCommSize;
	
	// Broadcasts the total size of R1, so that each process can allocate the full matrix to receive the data
	MPI_Bcast(&sizeR1, 1, MPI_INT, MASTER, MPI_COMM_WORLD);

	// Extracts the slice of R1 that belongs to the master process, so the Allgather operation can be performed
	rmatrix sliceR1(R1(Lb(R1, ROW), slice, Lb(R1, COL), Ub(R1, COL))), sliceR2;
	
	// Every process gathers the slices of R1	
	MPI_Allgather(R1, sliceR1, MPI_COMM_WORLD);

	sliceR1 = rmatrix();
	
	// Scatters R2 among the processes
	MPI_Scatter(R2, sliceR2, slice, rank, MASTER, MPI_COMM_WORLD);
	
	// Creates the matrices that will store the partial results
	rmatrix sliceD(Lb(sliceR2, ROW), Ub(sliceR2, ROW), Lb(sliceR2, COL), Ub(sliceR2, COL));
	rmatrix sliceR2temp(Lb(sliceR2, ROW), Ub(sliceR2, ROW), Lb(sliceR2, COL), Ub(sliceR2, COL));
	
	// Computes the partial results
	if(numThreads > 1){
		parcialD2R2temp(sliceR2, R1, sliceD, sliceR2temp, numThreads);
	} else {
		parcialD2R2temp(sliceR2, R1, sliceD, sliceR2temp);
	}
	
	// Gathers the partial results from the processes, including the master process
	MPI_Gather(D, sliceD, MASTER, MPI_COMM_WORLD);
	MPI_Gather(R2temp, sliceR2temp, MASTER, MPI_COMM_WORLD);

	return;
}

/**
 * calculateDR2temp - Function called by the SlaveLss class. All the slaves receive a slice of R2 matrix from the
 *                    master process. The part of R1 is the one that was received by the process when executing the
 *                    IminusAB function, so each process have to perform an MPI_Allgather operation to gather again
 *                    the R1 matrix.
 *            
 * Input parameters:
 *    rmatrix& sliceR1: slice of R1 matrix, reused from the IminusAB function;
 *    rmatrix& sliceR2: slice of R2 matrix that will be received from the master;
 *            int rank: process rank;
 *      int numThreads: number of threads that the computation should be subdivided into;
 * 
 * Output parameters:
 *  fstream& debugFile: reference to an already opened debug file.
 * 
 * Return value: void.
 * 
 */
void calculateDR2temp(rmatrix& sliceR1, rmatrix& sliceD, rmatrix& sliceR2temp, int rank, int numThreads){
	int sizeR1 = 0;
	
	// Receives the total size of R1
	MPI_Bcast(&sizeR1, 1, MPI_INT, MASTER, MPI_COMM_WORLD);
	
	rmatrix R1(sizeR1, sizeR1), sliceR2;
	
	// Gathers all the parts of R1 from the processes
	MPI_Allgather(R1, sliceR1, MPI_COMM_WORLD);
	// Receives the slice of R2
	MPI_Scatter(sliceR2, sliceR2, 0, rank, MASTER, MPI_COMM_WORLD);
	
	// Creates the matrices that will store the partial results
	sliceD = rmatrix(Lb(sliceR2, ROW), Ub(sliceR2, ROW), Lb(sliceR2, COL), Ub(sliceR2, COL));
	sliceR2temp = rmatrix(Lb(sliceR2, ROW), Ub(sliceR2, ROW), Lb(sliceR2, COL), Ub(sliceR2, COL));
	
	// Computes the partial results
	if(numThreads > 1){
		parcialD2R2temp(sliceR2, R1, sliceD, sliceR2temp, numThreads);
	} else {
		parcialD2R2temp(sliceR2, R1, sliceD, sliceR2temp);	
	}
	
	// Sends the partial results to the master process
	MPI_Gather(sliceD, sliceD, MASTER, MPI_COMM_WORLD);
	MPI_Gather(sliceR2temp, sliceR2temp, MASTER, MPI_COMM_WORLD);
	
	return;
}

/**
 * parcialD2R2temp - Function called by the calculateDR2temp. Performs the parcial computation.
 * 
 * Input parameters:
 *       rmatrix& a: left operand;
 *       rmatrix& b: right operand;
 * 
 * Output parameters:
 *       rmatrix& c: result of (a * b);
 *  rmatrix& r2temp: result of (a * b) - c, where (a * b) is computed using the long accumulator.
 * 
 * Return value: void.
 * 
 */
void parcialD2R2temp(rmatrix& a, rmatrix& b, rmatrix& c, rmatrix& r2temp){
	dotprecision ac(0);
	
	// The original version of LSS multiplies "a" and "b" using the "*" operator (e.g.: c = a * b).
	// Because multiplying matrices this way is not convenient to build a multi-threaded algorithm,
	// the code that performs this operation was extracted from the C-XSC source and the operation was
	// re-written
	for(int i = Lb(a, ROW); i <= Ub(a, ROW); i++){
		for(int j = Lb(b, COL); j <= Ub(b, COL); j++){
			ac = 0.0;
			accumulate(ac, a[i], b[Col(j)]);
			c[i][j] = rnd(ac);

			ac = ac - c[i][j];	
			r2temp[i][j] = rnd(ac);
		}
	}

	return;
}

/**
 * parcialD2R2temp - Function called by the calculateDR2temp.
 *                   Performs the parcial computation using nt threads.
 *            
 * Input parameters:
 *       rmatrix& a: left operand;
 *       rmatrix& b: right operand;
 *           int nt: number of threads;
 * 
 * Output parameters:
 *       rmatrix& c: result of (a * b);
 *  rmatrix& r2temp: result of (a * b) - c, where (a * b) is computed using the long accumulator.
 * 
 * Return value: void.
 * 
 */
void parcialD2R2temp(rmatrix& a, rmatrix& b, rmatrix& c, rmatrix& r2temp, int nt){
	int numThreads = nt, nRowsA = 0, threadSlice = 0;
	threadArgsDR2 args[numThreads];
	pthread_t threads[numThreads];
	pthread_attr_t attr;
	
	nRowsA = Ub(a, ROW) - Lb(a, ROW) + 1;
	threadSlice = nRowsA / nt;
	
	pthread_attr_init(&attr);
	pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
	
	for(int t = 0; t < numThreads; t++){
		args[t].threadId = t;
		args[t].sliceA = &a;
		args[t].b = &b;
		args[t].sliceC = &c;
		args[t].sliceR2temp = &r2temp;
		args[t].startRowA = Lb(a, ROW) + (t * threadSlice);
		args[t].endRowA = Lb(a, ROW) + ((t+1) * threadSlice) - 1;
		
		pthread_create(&threads[t], &attr, threadParcialDR2temp, &args[t]);
	}
	
	for(int t = 0; t < numThreads; t++){
		pthread_join(threads[t], NULL);
	}
	
	return;
}

/**
 * threadParcialDR2temp - Function executed by each thread.
 *            
 * Input parameters:
 *  void *tArgs: a pointer to the structure threadArgsDR2.
 * 
 * Output parameters: none.
 *  
 * Return value: void *
 * 
 */
void *threadParcialDR2temp(void *tArgs){
	threadArgsDR2 *args = (threadArgsDR2 *) tArgs;
	
	rmatrix &sliceA = *(args->sliceA);
	rmatrix &b = *(args->b);
	rmatrix &sliceC = *(args->sliceC);
	rmatrix &sliceR2temp = *(args->sliceR2temp);
	
	dotprecision ac(0);
	
	for(int i = args->startRowA; i <= args->endRowA; i++){
		for(int j = Lb(b, COL); j <= Ub(b, COL); j++){
			ac = 0.0;
			accumulate(ac, sliceA[i], b[Col(j)]);
			sliceC[i][j] = rnd(ac);
			
			ac = ac - sliceC[i][j];
			sliceR2temp[i][j] = rnd(ac);
		}
	}
	
	pthread_exit(NULL);
}
